On the Numerical Integration of Einstein's Field Equations 



OO 
On 



Thomas W. Baumgarte 1 and Stuart L. Shapiro 2 
1 Department of Physics, University of Illinois at Urbana- Champaign, Urbana, II 61801 
2 Departments of Physics and Astronomy and NCSA, University of Illinois at Urbana- Champaign, Urbana, II 61801 

Many numerical codes now under development to solve Einstein's equations of general relativity 
in 3+1 dimensional spacetimes employ the standard ADM form of the field equations. This form 
involves evolution equations for the raw spatial metric and extrinsic curvature tensors. Follow- 
ing Shibata and Nakamura, we modify these equations by factoring out the conformal factor and 
introducing three "connection functions". The evolution equations can then be reduced to wave 
equations for the conformal metric components, which are coupled to evolution equations for the 
connection functions. We evolve small amplitude gravitational waves and make a direct comparison 
of the numerical performance of the modified equations with the standard ADM equations. We find 
that the modified form exhibits much improved stability. 
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I. INTRODUCTION 

The physics of compact objects is entering a particu- 
larly exciting phase, as new instruments can now yield 
unprecedented observations. For example, there is evi- 
dence that the Rossi X-ray Timing Explorer has identi- 
fied the innermost stable circular orbit around an ac- 
creting neutron star Also, the new generation of 
gravitational wave detectors under construction, includ- 
ing LIGO, VIRGO, GEO and TAMA, promise to detect, 
for the first time, gravitational radiation directly (see, 
e.g., |). 

In order to learn from these observations (and, in the 
case of the gravitational wave detectors, to dramatically 
increase the likelihood of detection), one has to predict 
the observed signal from theoretical modeling. The most 
promising candidates for detection by the gravitational 
wave laser interferometers are the coalescences of black 
hole and neutron star binaries. Simulating such mergers 
requires self-consistent, numerical solutions to Einstein's 
field equations in 3 spatial dimensions, which is extremely 
challenging. While several groups, including two "Grand 
Challenge Alliances" S, have launched efforts to simu- 
late the coalescence of compact objects (see also 
the problem is far from being solved. 

Before Einstein's field equations can be solved nu- 
merically, they have to be cast into a suitable initial 
value form. Most commonly, this is done via the stan- 
dard 3+1 decomposition of Arnowitt, Deser and Misner 
(ADM, H). In this formulation, the gravitational fields 
are described in terms of spatial quantities (the spatial 
metric and the extrinsic curvature), which satisfy some 
initial constraints and can then be integrated forward in 
time. The resulting u g — K" equations are straightfor- 
ward, but do not satisfy any known hyperbolicity con- 
dition, which, as it has been argued, may cause stabil- 
ity problems in numerical implementations. Therefore, 
several alternative, hyperbolic formulations of Einstein's 
equations have been proposed [@-|l2| . Most of these for- 



mulations, however, also have disadvantages. Several of 
them introduce a large number of new, first order vari- 
ables, which take up large amounts of memory in numer- 
ical applications and require many additional equations. 
Some of these formulations require taking derivatives of 
the original equations, which may introduce further in- 
accuracies, in particular if matter sources are present. It 
has been widely debated if such hyperbolic formulations 
have computational advantages ]l3[ ]; their performance 
has yet to be compared directly with that of the origi- 
nal ADM equations. Accordingly, it is not yet clear if or 
how much the numerical behavior of the ADM equations 
suffers from their non-hyperbolicity. 

In this paper, we demonstrate by means of a numer- 
ical experiment and a direct comparison that the stan- 
dard implementation of the ADM system of equations, 
consisting of evolution equations for the bare metric and 
extrinsic curvature variables, is more susceptible to nu- 
merical instabilities than a modified form of the equa- 
tions based on a conformal decomposition as suggested 
by Shibata and Nakamura H] . We will refer to the stan- 
dard, "g — K" form of the equations as "System I" (see 
Section p A| below). We follow Shibata and Nakamura 
and modify these original ADM equations by factoring 
out a conformal factor and introducing a spatia l fiel d of 
connection functions ( "System II" , see Section II B be- 
low). The conformal decomposition separates "radiative" 
variables from "nonradiative" ones in the spirit of the 
"York-Lichnerowicz" split 15 18]. With the help of the 



connection functions, the Ricci tensor becomes an elliptic 
operator acting on the components of the conformal met- 
ric. The evolution equations can therefore be reduced to 
a set of wave equations for the conformal metric compo- 
nents, which are coupled to the evolution equations for 
the connection functions. These wave equations reflect 
the hyperbolic nature of general relativity, and can also 
be implemented numerically in a straight-forward and 
stable manner. 

We evolve low amplitude gravitational waves in pure 



1 



vacuum spacetimes, and directly compare Systems I and 
II for both geodesic slicing and harmonic slicing. We find 
that System II is not only more appealing mathemati- 
cally, but performs far better numerically than System 
I. In particular, we can evolve low amplitude waves in a 
stable fashion for hundreds of light travel timescales with 
System II, while the evolution crashes at an early time 
in System I, independent of gauge choice. We present 
these results in part to alert developers of 3+1 general 
relativity codes, many of whom currently employ System 
I, that a better set of equations may exist for numerical 
implementation. 

The paper is organized as follows. In Section [H|, we 
present the basic equations of both Systems I and II. 
We briefly discuss our numerical implementation in Sec- 
tion and present numerical results in Section IV . In 
Section 0, we summarize and discuss some of the impli- 
cations of our findings. 



Here Di is the covariant derivative associated with 7^, 
Rij is the three-dimensional Ricci tensor 

R*3 — ~jl kl {lkjM + lil,kj ~ lkl,ij ~ lij,kl) CO 

+i kl ( r ™r mfei - T™T mk i^ , 

and R is its trace R = j^Rij. We have also introduced 
the matter sources p, Si and Sij , which are projections of 
the stress-energy tensor with respect to the unit normal 
vector n a , 



p = n a nf3T a(3 , 
Si = - lia n p T a \ 
Sij = ^i a ljf}T aP , 
and have abbreviated 



(8) 



II. BASIC EQUATIONS 

A. System I 

We write the metric in the form 

ds 2 = -u 2 dt 2 + "f tj (dx l + ^dt){dx j + ftdt), (1) 

where a is the lapse function, (3 1 is the shift vector, and 
7y is the spatial metric. Throughout this paper, Latin 
indices are spatial indices and run from 1 to 3, whereas 
Greek indices are spacetime indices and run from to 
3. The extrinsic curvature can be defined by the 
equation 



where 



— jij = -2aKij, 



dt ~ dt 



(2) 



(3) 



and where Cp denotes the Lie derivative with respect to 

The Einstein equations can then be split into the 
Hamiltonian constraint 



R - KijK ij + K 2 = 2p, 
the momentum constraint 



(4) 



D 3 K\ - DiK = Si, (5) 
and the evolution equation for the extrinsic curvature 
d 



dt 



Kij = -DiDjOt + a(R 2J - 2K a K l j + KK 2j - M v 



(0) 



M i:j = S i:j + -lij(p-S), 



(9) 



where S is the trace of Sij , S = 7 y Sij . 

The evolution equations (||) and (||) together with the 
constraint equations (0) and (||) are equivalent to the 
Einstein equations, and are commonly referred to as the 
ADM form of the gravitational field equations ||[l7| . We 
will call these equations System I. This system is widely 
used in numerical relativity calculations (e.g. p8| , p9| ), 
even though its mathematical structure is not simple to 
characterize and may not be ideal for computation. In 
particular, the Ricci tensor (f7j) is not an elliptic opera- 
tor: while the last one of the four terms involving second 
derivatives, j kl ~fij,ki, is an elliptic operator acting on the 
components of the metric, the elliptic nature of the whole 
operator is spoiled by the other three terms involving sec- 
ond derivatives. Accordingly, the system as a whole does 
not satify any known hyperbolicity condition (see also 
the discussion in Jll|). Therefore, to establish existence 
and uniqueness of solutions to Einstein's equations, most 
mathematical analyses rely either on particular coordi- 
nate choices or on different formulations. 



B. System II 



Instead of evolving the metric 7^ and the extrinsic 
curvature K^, we can evolve a conformal factor and 
the trace of the extrinsic curvature separately ("York- 
Lichnerowicz split" [ ^5|Jl^ ]). Such a split is very appeal- 
ing from both a theoretical and computational point of 
view, and has been widely applied in numerical axisym- 
metric (2+1) calculations (sec, e.g., Q). More recently, 
Shibata and Nakamura J14| applied a similar technique 
in a three-dimensional (3+1) calculation. Adopting their 
notation, we write the conformal metric as 



(10) 
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and choose 



e 4 * = 7 1/8 = det( 7 «) 1/3 . 



(11) 



so that the determinant of %j is unity. We also write the 
trace-free part of the extrinsic curvature K%j as 



Aij — Kij ^jijK, 



(12) 



where K = j^Kij. It turns out to be convenient to 
introduce 



At 



(13) 



where Di is the derivative operator associated with 7^ , 
and D l = fWj. 

The "tilde" Ricci tensor i?y is the Ricci tensor asso- 
ciated with 7^ and could be computed by inserting 7^ 
into equation (pi). However, we can bring the Ricci tensor 
into a manifestly elliptic form by introducing the "con- 
formal connection functions" 



~ l 1 

■7 . 



(21) 



where the T^ k are the connection coefficients associated 
with jij, and where the last equality holds because 7 = 1. 
In terms of these, the Ricci tensor can be written [p22| 



We will raise and lower indices of A^ with the conformal 
metric %, so that A* = e^A^ (see @). 

Taking the trace of the evolution equations (|J) and (|^) 
with respect to the physical metric 7y , we find [M 



— d> — — aK 
dV 6 



(14) 



and 
d 



-K 



rJi , . -flDjDia + afaA* + ^K 2 ) + ±a(p + S), 

(15) 



where we have used the Hamiltonian constraint (Q) to 
eliminate the Ricci scalar from the last equation. The 
tracefree parts of the two evolution equations yield 



dt^ 



-2aAij. 



(16) 



and 
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(D iDj af F + a{Rj 3 F - fig*)) 



-a{KA lJ -2A il A l j ). 



(17) 



In the last equation, the superscript TF denotes the 
trace- free part of a tensor, e.g. Rjj F — Rij — jijR/3. 
Note that the trace R could again be eliminated with 
the Hamiltonian constraint (Q). Note also that 7y and 
Aij are tensor densities of weight —2/3, so that their Lie 
derivative is, for example, 

CpAii = /? fc d fe iy + A lk d 3 (3 k + A kj dif3 k - ^A l0 d k f3 k . 

(18) 

The Ricci tensor R^ in ( |i"7| ) can be written as the sum 



Rij - !>'., I />',,■ 



(19) 



Here Rf- is 



R% = -2DiDj^> - 2%D l D i( j> 

+i{b i 4,){b j 4,) - 4 7y pV)(A0), (20) 



Rij = --jl lm lij,lm + lk(id 3 )t k + f fe f( ii)fe + 

l lm ( 2 f f(if j) km + f k m T k ij 



(22) 



The principal part of this operator, j lm jij,i m , is that of a 
Laplace operator acting on the components of the metric 
7y. It is obviously elliptic and diagonally dominant (as 
long as the metric is diagonally dominant). All the other 
second derivatives of the metric appearing in (Q) have 
been absorbed in the derivatives of the connection func- 
tions. At least in appropriately chosen coordinate sys- 
tems (for example j3 l = 0), equations |jl]) and (Jl?]) there- 
fore reduce to a coupled set of nonlinear, inhomogeneous 
wave equations for the conformal metric 7y , in which the 
gauge terms K and T l , the conformal factor exp ((/>), and 
the matter terms Mij appear as sources. Wave equa- 
tions not only reflect the hyperbolic nature of general 
relativity, but can also be implemented numerically in a 
straight-forward and stable manner. The same method 
has often been used to reduce the four-dimensional Ricci 
tensor R a p [ p3[ and to brin g E instein's equations into a 
symmetric hyperbolic form p4[ . 

Note that the connection functions L l are pure gauge 
quantities in the sense that they could be chosen, for 
example, to vanish by a suitable choice of spatial coor- 
dinates ("conformal three- harmonic coordinates", com- 
pare 1 25 1). The P would then play the role of "con- 
formal gauge source functions" (compare |2^,Q). Here, 
however, we impose the gauge by choosing the shift 
and evolve the T l with equation ( |24| ) below. Similarly, K 
is a pure gauge variable (and could be chosen to vanish 
by imposing maximal time slicing). 

An evolution equation for the T l can be derived by 
permuting a time derivative with the space derivative 
in © 



dt 



+ 37V 



/?<7 



(23) 



It turns out to be essential for the numerical stability of 
the system to eliminate the divergence of A 1 ^ with the 
help of the momentum constraint (pf), which yields 
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J^F = -2A ij a tj + 2a{f) k A kj - 



JL (pl^^ _ 2^(3%, + \fip 1 ^ . (24) 



We now consider <j), K, 7y, Aij and T 1 as fundamen- 
tal variables. These can be evolved with the evolution 
equations @, (|l|), |L6|), (|l7|), and (|§), which we call 



System II. Note that obviouly not all these variables are 
independent. In particular, the determinant of 7^ has 
to be unity, and the trace of Aij has to vanish. These 
conditions can either be used to reduce the number of 
evolved quantities, or, alternatively, all quantities can be 
evolved and the conditions can be used as a numerical 
check (which is what we do in our implementation). 



III. NUMERICAL IMPLEMENTATION 

In order to compare the properties of Systems I and 
II, we implemented them numerically in an identical en- 
vironment. We integrate the evolution equations with 
a two-level, iterative Crank-Nicholson method. The it- 
eration is truncated after a certain accuracy has been 
achieved. However, we iterate at least twice, so that the 
scheme is second order accurate. 

The gridpoints on the outer boundaries are updated 
with a Sommerfeld condition. We assume that, on the 
outer boundaries, the fundamental variables behave like 
outgoing, radial waves 



Q(t,r) 



G{at - e 2rt >r) 



(25) 



Here Q is any of the fundamental variables (except for 
the diagonal components of 7y , for which the radiative 
part is Q — ju — 1), and G can be found by following 
the characteristic back to the previous timestep and in- 
terpolating the corresponding variable to that point (see 
also [Q). We found that a linear interpolation is ade- 
quate for our purposes. 

We impose octant symmetry in order to minimize the 
number of gridpoints, and impose corresponding symme- 
try boundary conditions on the symmetry plains. Unless 
noted otherwise, the calculations presented in this paper 
were performed on grids of (32) 3 gridpoints, and used a 
Courant factor of 1/4. The code has been implemented 
in a parallel environment on SGI Power ChallengeArray 
and SGI CRAY Origin2000 computer systems at NCSA 
using DAGH p(| software for parallel processing. 



IV. RESULTS 
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FIG. 1. Evolution of the trace of the extrinsic curvature K 
for a small amplitude wave in geodesic slicing at the origin 
(see text for details). The solid line is the result for System 
II, and the dashed line for System I. The dotted line is the 
approximate solution ([291). 



A. Initial Data 

For initial data, we choose a linearized wave solution 
(which is then evolved with the full nonlinear systems I 
and II). Following Teukolsky ^7j, we construct a time- 
symmetric, even-parity L = 2, M = solution. The 
coefficients A, B and C (see equation (6) in |27]]) are de- 
rived from a function 



F(t, r)=A(t±r) exp(-A(< ± rf 



(26) 



Unless noted otherwise, we present results for an am- 
plitude A = 10~ 3 and a wavelength A = 1. The outer 
boundary conditions are imposed at x,y,z = 4. 
We evolve these initial data for zero shift 



l3 l = 0, 



(27) 



and compare the performance of Systems I and II for 
both geodesic and harmonic slicing. 



B. Geodesic Slicing 



In geodesic slicing, the lapse is unity 
a = 1. 



(28) 



Since the acceleration of normal observers satisfies a a — 
D a In a = 0, these observers follow geodesies. The energy 
content of even a small, linear wave packet will therefore 
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FIG. 2. Evolution of the extrinsic curvature component 
K zz at the origin in geodesic slicing. The solid line is the 
result for System II, and the dashed line for System I. For 



System II, we constructed K zz from A z 



K and j z 



focus these observers, and even after the wave has dis- 
persed, the observers will continue to coast towards each 
other. Since (3 l = 0, normal observers are identical to co- 
ordinate observers, hence geodesic slicing will ultimately 
lead to the formation of a coordinate singularity even for 
arbitrarily small waves. 

The timescale for the formation of this singularity can 
be estimated from equation ( |l5| ) with a = 1 and p l — 0. 
The Aij, which can be associated with the gravitational 
waves, will cause K to increase to some finite value, say 
Kq at time to, even if K was zero initially. After roughly 
a light crossing time, the waves will have dispersed, and 
the further evolution of K is described by dtK ~ K 2 /3, 



K 



3K 



K (t - t a ) 



(29) 



(see |14)). Obviously, the coordinate singularity forms at 
t ~ 3 /Kq + to as a result of the nonlinear evolution. 

We can now evolve the wave initial data with Systems 
I and II and compare how well they reproduce the for- 
mation of the coordinate singularity. 

In Figure 1, we show K at the origin (x = y = z = 0) 
as a function of time both for System I (dashed line) 
and System II (solid line). We also plot the approx- 
imate analytic solution (B9) as a dotted line, which 
we have matched to the System I solution with values 
Kq = 0.00518 and t = 10. For these values, equa- 
tion ( p9|) predicts that the coordinate singularity appears 
at t ~ 590. In the insert, we show a blow-up of System II 
for early times. It can be seen very clearly how the initial 



wave content lets K grow from zero to the "seed" value 
Ko- Once the waves have dispersed, System II approx- 
imately follows the solution ( p^ ) up to fairly late times. 
System I, on the other hand, crashes long before the co- 
ordinate singularity appears. 

In Figure 2, we compare the extrinsic curvature com- 
ponent K zz evaluated at the origin. The noise around 
t ~ 8, which is present in the evolutions of both systems, 
is caused by reflections of the initial wave off the outer 
boundaries. It is obvious from these plots that System 
II evolves the equations stably to a fairly late time, at 
which the integration eventually becomes inaccurate as 
the coordinate singularity approaches. We stopped this 
calculation when the iterative Crank-Nicholson scheme 
no longer converged after a certain maximum number of 
iterations. It is also obvious that System I performs ex- 
tremely poorly, and crashes at a very early time, well 
before the coordinate singularity. 

It is important to realize that the poor performance of 
System I is not an artifact of our numerical implementa- 
tion. For example, the ADM code currently being used 
by the Black Hole Grand Challenge Alliance, is based on 
the equations of System I, and also crashes after a very 
similar time (see also |Q, where a run with a much 
smaller initial amplitude nevertheless crashes earlier than 
our System II). This shows that the code's crashing is 
intrinsic to the equations and slicing, and not to our nu- 
merical implementation. 



C. Harmonic Slicing 

Since geodesic slicing is known to develop coordinate 
singularities for generic, nontrivial initial data, it is ob- 
viously not a very good slicing condition. We therefore 
also compare the two Systems using harmonic slicing. 
In harmonic slicing, the coordinate time t is a harmonic 
function of the coordinates V Q V Q i = 0, which is equiva- 
lent to the condition 



r° 



9 1 a/3 







(30) 



(where the are the connection coefficients associated 
with the four-dimensional metric g a p). For f3 l = 0, the 
above condition reduces to 

d t a = -a 2 K. (31) 

Inserting (|l4|), this can be written as 

<9 t (ae- 60 )=O or a = C(^)e 6 *, (32) 

where C'(x l ) is a constant of integration, which depends 
on the spatial coordinates only. In practice, we choose 
C{x l ) = 1. 

In Figure 3, we show results for the same initial data 
as in the last section. Obviously, both Systems do much 
better for this slicing condition. System I crashes much 
later than in geodesic slicing (after about 40 light crossing 
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times. While efforts have been undertaken to stabilize 
such codes with the help of appropriate outer boundary 
conditions fl8| , p9| ], our findings point to the equations 
themselves as the fundamental cause of the problem, and 
not to the outer boundaries. Obviously, boundary condi- 
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t 

FIG. 3. Evolution of the extrinsic curvature component 
K zz at the origin in harmonic slicing. The solid line is the 
result for System II, and the dashed line for System I. For 
System II, we constructed K zz from A zz , 4>, K and 7 ZZ . 



times, as opposed to about 10 for geodesic slicing), but 
it still crashes. System II, on the other hand, did not 
crash after even over 100 light crossing times. We never 
encountered a growing instability that caused the code 
to crash. 



V. SUMMARY AND CONCLUSION 

We numerically implement two different formulations 
of Einstein's field equations and compare their perfor- 
mance for the evolution of linear wave initial data. Sys- 
tem I is the standard set of ADM equations for the evolu- 
tion of 7y- and . In System II, we conformally decom- 
pose the equations and introduce connection functions. 
The conformal decomposition naturally splits "radiative" 
variables from "nonradiative" ones, and the connection 
functions are used to bring the Ricci tensor into an ellip- 
tic form. These changes are appealing mathematically, 
but also have a striking numerical consequence: System 
II performs far better than System I. 

It is interesting to note that most earlier axisymmetric 
codes (e.g. [^(J) also relied on a decomposition similar to 
that of System II. Much care was taken to identify ra- 
diative variables and to integrate those variables as op- 
posed to the raw metric components. It is surprising 
that this experience was abandoned in the development 
of most 3+1 codes, which integrate equations equivalent 
to System I. These codes have been partly successful [[l9| , 
but obvious problems remain, as for example the inabil- 
ity to integrate low amplitude waves for arbitrarily long 



tions as employed in the perturbative approach in [|18 29 
or in the characteristic approach in p0| are still needed 
for accuracy - but our results clearly suggest that they 
are not needed for stability ]3l| ]. 

Some of the recently proposed hyperbolic systems are 
very appealing in that they bring the equations in a first 
order, symmetric hyperbolic form, and that all character- 
istics are physical (i.e., are either at rest with respect to 
normal observers or travel with the speed of light) [p|,^2[. 
These properties may be very advantageous for numerical 
implementations, in particular at the boundaries (both 
outer boundaries and, in the case of black hole evolutions, 
inner "apparent horizon" boundaries). Some of these sys- 
tems have also been implemented numerically, and show 
stability properties very similar to our System II p2| . 
Our System II, on the other hand, uses fewer variables 
than most of the hyperbolic formulations, and does not 
take derivatives of the equations, which may be advanta- 
geous especially when matter sources are present. This 
suggests that a system similar to System II may be a good 
choice for evolving interior solutions and matter sources, 
while one may want to match to one of the hyperbolic 
formulations for a better treatment of the boundaries. 

The mathematical structure of System II is more ap- 
pealing than that of System I, and these improvements 
are reflected in the numerical behavior. We therefore con- 
clude that the mathematical structure has a very deep 
impact on the numerical behavior, and that the ability 
to finite difference the standard u g — K" ADM equations 
may not be sufficient to warrant a stable evolution. 
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